/**
 * @file mt2d_edge.cpp
 * @author Yi Zhang (zhangyiss@icloud.com)
 * @brief 2-D Forward modeling of magnetotelluric data at a given frequency using an edge_based finite-element method
 * @version 1.0
 * @date 2021-12-22
 * 
 * @copyright Copyright (c) 2021
 * 
 * 
 * +z downward, z=0 the Air-Earth interface
 * time dependence exp(-iwt)
 * 
 * TE-mode: E = (Ex, 0, 0), and we solve H = (0, Hy, Hz)
 * TM-mode: H = (Hx, 0, 0), and we solve E = (0, Ey, Ez)
 *          x
 *        /
 *      /
 *    /    air (conductivity = 1e-10, epsilon = 1.0, mu = 1.0, bottom depth = 0)
 *  /
 * O---------------------> y
 * |   ...
 * |---------------------
 * |   layer bottom (conductivity = 0.1, epsilon = 1.0, mu = 1.0, bottom depth = 1e+20)
 * |
 * z
 * 
 */

/*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
#if __INTELLISENSE__
#undef __ARM_NEON
#undef __ARM_NEON__
#endif

// Provide FEM shape functions and the Gauss Legendre Numerical Integral (GLNI) functions
#include "../utilities/glni.h"
#include "../utilities/efem_sf.h"
// Provide primary filed and boundary values
#include "../utilities/layerEarth_MT.h"
// Provide array structures and file io
#include "gctl/core.h"
#include "gctl/io.h"
// Provide sparse matrix and solver
#include "Eigen/Dense"
#include "Eigen/Sparse"

#define pi 3.1415926535897932384626433832795
#define _mu0 (4e-7*pi)
#define _eps0 8.8541878176e-12

typedef std::complex<double> complex_d;
typedef Eigen::SparseMatrix<complex_d> spmat_cd; // 注意Eigen默认的稀疏矩阵排序为列优先
typedef Eigen::Triplet<complex_d> triplt_cd;
typedef Eigen::VectorXcd vector_cd;

static const complex_d _one = complex_d(1.0, 0.0);
static const complex_d _zero = complex_d(0.0, 0.0);
static const complex_d _1i  = complex_d(0.0, 1.0);

enum MT_Mode {TE, TM, Modeling};

struct conductive_region
{
    int tag;
    double cndt, mu, epsl; // cndt (conductivity) epsl (epsilon)
    double ref_cndt; // background conductivity
};

struct mt_station
{
    gctl::vertex2dc *vert;
    double app_resist; // apparent resistivity
    double phase; // phase angle
    std::vector<gctl::triangle2d*> host_tris;
};

struct edge_setup
{
    unsigned int uv_id[2];
    efem_sf::ordertype_e uv_type[2];
};

class mt2d
{
public:
    mt2d(){}
    virtual ~mt2d(){}
    void run(int argc, char *argv[]); // 执行计算流程

protected:
    void read_config(std::string filename); // 读入配置文件
    void read_model(std::string filename); //读入Triangle模型
    void write_model(std::string filename, bool debug = true); // 写入Gmsh模型文件
    void write_observation(std::string filename);

    double get_element_conductivity(int tag_num, bool ref_cndt = false);
    void assemble_linear_system();
    void apply_boundary_condition();
    void solve();
    void post_process();

    static double _Aij_part1(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance)
	{
		return reinterpret_cast<mt2d*>(instance)->Aij_part1(x, y, x1, x2, x3, y1, y2, y3, att);
	}

	double Aij_part1(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att) // 偏微分方程的左端项1待积分函数
	{
		edge_setup *ft = static_cast<edge_setup*>(att);
        double curl_u, curl_v, dummy;

        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[0], efem_sf::Curl, ft->uv_type[0], curl_u, dummy);
        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[1], efem_sf::Curl, ft->uv_type[1], curl_v, dummy);
        return curl_u*curl_v;
	}

    static double _Aij_part2(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance)
	{
		return reinterpret_cast<mt2d*>(instance)->Aij_part2(x, y, x1, x2, x3, y1, y2, y3, att);
	}

    double Aij_part2(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att) // 偏微分方程的左端项2待积分函数
	{
		edge_setup *ft = static_cast<edge_setup*>(att);
        double ux, uy, vx, vy;

        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[0], efem_sf::Value, ft->uv_type[0], ux, uy);
        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[1], efem_sf::Value, ft->uv_type[1], vx, vy);
        return (ux*vx + uy*vy);
	}

    // 偏微分方程的右端项
    static std::complex<double> _Bi_TM(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance)
	{
		return reinterpret_cast<mt2d*>(instance)->Bi_TM(x, y, x1, x2, x3, y1, y2, y3, att);
	}
    
    std::complex<double> Bi_TM(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att) // 偏微分方程的右端项待积分函数
	{
        edge_setup *ft = static_cast<edge_setup*>(att);
        double ux, uy;

        le_.calculate_EH(-1.0*y, x);

        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[0], efem_sf::Value, ft->uv_type[0], ux, uy);

        // 这里原本应该采用一次场在棱边上的分量进行一个双循环积分，即 Ni*Nj*Ej
        // 因为 \sum_e{Nj*Ej} 等于一次场在积分点处的场值 因此这里转换成了使用一次场值的单循环积分
        // 参考文献
        // Cai, H., Xiong, B., Han, M., & Zhdanov, M. (2014). 3D controlled-source electromagnetic modeling 
        // in anisotropic medium using edge-based finite element method. Computers & Geosciences, 73, 164–176. 
        // https://doi.org/10.1016/j.cageo.2014.09.008
        return (ux*le_.get_Ey() - uy*le_.get_Ez());
	}

    static std::complex<double> _Bi_TE(double x, double y, double x1, double x2, double x3, 
        double y1, double y2, double y3, void *att, void *instance)
    {
        return reinterpret_cast<mt2d*>(instance)->Bi_TE(x, y, x1, x2, x3, y1, y2, y3, att);
    }
    
    std::complex<double> Bi_TE(double x, double y, double x1, double x2, double x3, 
        double y1, double y2, double y3, void *att) // 偏微分方程的右端项待积分函数
    {
        edge_setup *ft = static_cast<edge_setup*>(att);
        double curl_u, dummy;

        le_.calculate_EH(-1.0*y, x);

        elsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, ft->uv_id[0], efem_sf::Curl, ft->uv_type[0], curl_u, dummy);
        
        // 这里将右端项视为由一次场与异常电导率共同构成的场源项 注意y轴（竖轴）的方向与层状模型的z轴方向相反 
        // 参考文献
        // Key, K., & Weiss, C. (2006). Adaptive finite-element modeling using unstructured grids: 
        // The 2D magnetotelluric example. Geophysics, 71(6), G291–G299. https://doi.org/10.1190/1.2348091
        return -1.0*curl_u*le_.get_Ex();
    }

    // 以下为类的变量
    unsigned int elem_num_; // 有限元个数（分割的单元数）
    unsigned int node_num_; // 节点个数
    unsigned int edge_num_;
    gctl::array<int> node_markers_; // 边界顶点的标签值
    gctl::array<int> elem_markers_; // 元素的标签值
    gctl::array<gctl::vertex2dc> mesh_nodes_; // 网格顶点
    gctl::array<gctl::triangle2d> mesh_elems_; // 网格元素
    gctl::array<gctl::edge2d> mesh_edges_;
    gctl::matrix<gctl::edge2d*> elem_edge_ptr_;
    gctl::gmshio meshio_; // 网格文件IO

    glni::triangle g2d_tri_; // 二维三角形区域高斯勒让德积分对象
    efem_sf::linear elsf_; // 线性边矢量形函数对象
    spmat_cd A;
    vector_cd X, B, P;

    std::string tri_file_; // 网格文件名
    double omega_; // 电磁波角频率
    double theta_in_; // 主场入射角
    int station_marker_; // 观测点标签值
    int bound_marker_; // 边界顶点标签值 (我们并不需要这个参数，保留它只是为了和基于顶点的有限元方法共用一套参数文件)
    int air_marker_; // 空气层标签值
    std::vector<conductive_layer> primary_model_; // 主场模型 (导电率 相对磁导率 相对介电常数 底深度)
    gctl::array<conductive_region> secondary_model_; // 磁场模型 (标签值 导电率 相对磁导率 相对介电常数)

    int obs_num_;
    std::vector<mt_station> obs_site_; // 观测站

    layerEarthMT le_;

    MT_Mode work_mode_;
};

void mt2d::read_config(std::string filename)
{
    std::clog << "Reading configuration files ... ";
    std::ifstream infile;
    gctl::open_infile(infile, filename, ".config");

    double frequency;
    std::string w_mode;
    infile >> tri_file_ >> w_mode >> frequency >> theta_in_ 
        >> station_marker_ >> bound_marker_ >> air_marker_;

    if (w_mode == "TE") work_mode_ = TE;
    else if (w_mode == "TM") work_mode_ = TM;
    else if (w_mode == "Modeling") work_mode_ = Modeling;
    else throw gctl::runtime_error("Invalid MT surveying mode.");

    omega_ = 2.0*GCTL_Pi*frequency;

    int tmp_size;
    infile >> tmp_size;
    if (tmp_size < 2)
    {
        throw gctl::runtime_error("Invalid primary model sizes: "+std::to_string(tmp_size));
    }

    // Set primary model
    primary_model_.resize(tmp_size);
    for (int i = 0; i < tmp_size; ++i)
    {
        infile >> primary_model_[i].cndt >> primary_model_[i].mu 
            >> primary_model_[i].epsl >> primary_model_[i].btm_depth;
    }

    le_.set_layers(primary_model_);
    if (work_mode_ == TM) le_.set_to_TM_mode();
    le_.set_frequency(frequency);
    le_.initialize();

    infile >> tmp_size;
    if (tmp_size <= 0)
    {
        throw gctl::runtime_error("Invalid secondary model sizes: "+std::to_string(tmp_size));
    }

    secondary_model_.resize(tmp_size);
    for (int i = 0; i < tmp_size; ++i)
    {
        infile >> secondary_model_[i].tag >> secondary_model_[i].cndt >> secondary_model_[i].ref_cndt 
            >> secondary_model_[i].mu >> secondary_model_[i].epsl;
    }

    infile.close();
    std::clog << "done\n";
    return;
}

// Read Triangle's output files
void mt2d::read_model(std::string filename)
{
    std::clog << "Reading Triangle files ... ";
    gctl::matrix<double> tmp_markers;
    gctl::read_Triangle_node(filename, mesh_nodes_, gctl::Packed, &node_markers_);
    gctl::read_Triangle_element(filename, mesh_elems_, mesh_nodes_, gctl::Packed, &tmp_markers);
    gctl::read_Triangle_neighbor(filename, mesh_elems_, gctl::Packed);

    elem_num_ = mesh_elems_.size();
    node_num_ = mesh_nodes_.size();

    elem_markers_.resize(elem_num_);
    for (int i = 0; i < elem_num_; ++i)
    {
        elem_markers_[i] = (int) tmp_markers[i][0];
    }
    tmp_markers.clear();

    // find edge neighbors
    std::vector<std::vector<int>> edge_n;
    edge_n.resize(elem_num_);
    for (int i = 0; i < elem_num_; i++)
    {
        edge_n[i].reserve(6); // 每个边的list最多6个
    }

    int small_id, big_id;
    for (int i = 0; i < elem_num_; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            if (mesh_elems_[i].neigh[j] != nullptr) // 邻居存在 公共边
            {
                // 这里我们总是按小-大排序序号，保证不会输入混乱
                small_id = std::min(mesh_elems_[i].id, mesh_elems_[i].neigh[j]->id);
                big_id   = std::max(mesh_elems_[i].id, mesh_elems_[i].neigh[j]->id);
                edge_n[small_id].push_back(big_id);
            }
            else // 邻居不存在 边界边
            {
                small_id = mesh_elems_[i].id;
                // 注意这里我们使用了一个小技巧来识别边界上的边，负值对应的边即为边界边
                // 为了避开0所以加了1
                big_id = -1*(j+1);
                edge_n[small_id].push_back(big_id);
            }
        }
    }

    // 去除edge_n中的重复输入
    std::vector<int>::iterator pos;
    for (int i = 0; i < elem_num_; i++)
    {
        if (!edge_n[i].empty()) // 很有可能某个边没有输入
        {
            std::sort(edge_n[i].begin(), edge_n[i].end()); //对顶点序列由小到大排序
            pos = std::unique(edge_n[i].begin(), edge_n[i].end()); //获取重复序列开始的位置
            edge_n[i].erase(pos, edge_n[i].end()); //删除重复点
        }
    }

    // 保存三角形的边的列表 首先统计边的数量
    edge_num_ = 0;
    for (int i = 0; i < elem_num_; i++)
        edge_num_ += edge_n[i].size();
    mesh_edges_.resize(edge_num_);

    int count = 0;
    gctl::vertex2dc *tmp_vert;
    for (int i = 0; i < elem_num_; i++)
    {
        if (!edge_n[i].empty())
        {
            for (int j = 0; j < edge_n[i].size(); j++)
            {
                mesh_edges_[count].id = count;
                mesh_edges_[count].neigh[0] = mesh_elems_.get(i);
                if (edge_n[i][j] >= 0) // 公共边 获取顶点索引并将邻居设置为对应指针
                {
                    mesh_edges_[count].neigh[1] = mesh_elems_.get(edge_n[i][j]);
                    // 如果相邻的两个三角形中的其中一个三角形的一个顶点不等于另一个三角形的所有顶点，则这个三角形的
                    // 另两个顶点连成的边即为公共边
                    for (int k = 0; k < 3; k++)
                    {
                        tmp_vert = mesh_elems_[i].vert[k];
                        if (tmp_vert != mesh_edges_[count].neigh[1]->vert[0] &&
                            tmp_vert != mesh_edges_[count].neigh[1]->vert[1] &&
                            tmp_vert != mesh_edges_[count].neigh[1]->vert[2])
                        {
                            mesh_edges_[count].vert[0] = mesh_elems_[i].vert[(k+1)%3];
                            mesh_edges_[count].vert[1] = mesh_elems_[i].vert[(k+2)%3];
                            break;
                        }
                    }
                }
                else // 边界边 获取顶点索引并将邻居设置为空
                {
                    mesh_edges_[count].neigh[1] = nullptr;
                    mesh_edges_[count].vert[0] = mesh_elems_[i].vert[(-1*edge_n[i][j])%3]; // 定位边的顶点索引
                    mesh_edges_[count].vert[1] = mesh_elems_[i].vert[(-1*edge_n[i][j]+1)%3];
                }
                count++;
            }
        }
    }

    //清除向量
    gctl::destroy_vector(edge_n);

    //搜索三角形各个边的全局索引
    elem_edge_ptr_.resize(elem_num_, 3, nullptr);
    for (int i = 0; i < edge_num_; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            if (mesh_edges_[i].neigh[j] != nullptr)
            {
                count = mesh_edges_[i].neigh[j]->id;
                for (int k = 0; k < 3; k++)
                {
                    if ((mesh_elems_[count].vert[k] == mesh_edges_[i].vert[0] &&
                        mesh_elems_[count].vert[(k+1)%3] == mesh_edges_[i].vert[1]) || 
                        (mesh_elems_[count].vert[k] == mesh_edges_[i].vert[1] &&
                        mesh_elems_[count].vert[(k+1)%3] == mesh_edges_[i].vert[0])) // 定位边在三角形中的局部索引
                    {
                        elem_edge_ptr_[count][k] = mesh_edges_.get(i);
                        break;
                    }
                }
            }
        }
    }

    std::clog << "done\n";
    std::clog << "triangle number = " << elem_num_ << ", node number = " << node_num_ << ", edge number = " << edge_num_ << std::endl;
    return;
}

// Output a Gmsh file
void mt2d::write_model(std::string filename, bool debug)
{
    std::clog << "Writing Gmsh file ... ";
    meshio_.init_file(filename, gctl::Output);
    meshio_.set_packed(gctl::NotPacked, gctl::Output);
    meshio_.save_mesh(mesh_elems_, mesh_nodes_);
    meshio_.save_data("vertex marker", node_markers_, gctl::NodeData);
    meshio_.save_data("element marker", elem_markers_, gctl::ElemData);

    if (debug)
    {
        gctl::array<double> total_area(node_num_, 0.0);
        gctl::array<double> d1(node_num_, 0.0), d2(node_num_, 0.0);
        gctl::array<double> d3(node_num_, 0.0), d4(node_num_, 0.0);
        double ux, uy;
        complex_d fx, fy;
        efem_sf::ordertype_e o_type;

        double tmp_area;
        for (size_t e = 0; e < elem_num_; e++)
        {
            // 双重加权平均 同时考虑三角形面积与电导率的权重 减少空气层界面的畸变
            tmp_area = mesh_elems_[e].area() * get_element_conductivity(elem_markers_[mesh_elems_[e].id]);

            for (size_t v = 0; v < 3; v++)
            {
                total_area[mesh_elems_[e].vert[v]->id] += tmp_area;
            }
        }

        for (size_t e = 0; e < elem_num_; e++)
        {
            tmp_area = mesh_elems_[e].area() * get_element_conductivity(elem_markers_[mesh_elems_[e].id]);

            for (size_t v = 0; v < 3; v++)
            {
                fx = fy = _zero;
                for (size_t d = 0; d < 3; d++)
                {
                    if (mesh_elems_[e].vert[d]->id < mesh_elems_[e].vert[(d+1)%3]->id) o_type = efem_sf::Normal;
                    else o_type = efem_sf::Reverse;

                    elsf_.triangle(mesh_elems_[e].vert[v]->x, mesh_elems_[e].vert[v]->y, 
                        mesh_elems_[e].vert[0]->x, mesh_elems_[e].vert[1]->x, mesh_elems_[e].vert[2]->x, 
                        mesh_elems_[e].vert[0]->y, mesh_elems_[e].vert[1]->y, mesh_elems_[e].vert[2]->y, 
                        d, efem_sf::Value, o_type, ux, uy);

                    fx += ux*X[elem_edge_ptr_[e][d]->id];
                    fy += uy*X[elem_edge_ptr_[e][d]->id];
                }

                d1[mesh_elems_[e].vert[v]->id] += fx.real()*tmp_area/total_area[mesh_elems_[e].vert[v]->id];
                d2[mesh_elems_[e].vert[v]->id] += fx.imag()*tmp_area/total_area[mesh_elems_[e].vert[v]->id];
                d3[mesh_elems_[e].vert[v]->id] += fy.real()*tmp_area/total_area[mesh_elems_[e].vert[v]->id];
                d4[mesh_elems_[e].vert[v]->id] += fy.imag()*tmp_area/total_area[mesh_elems_[e].vert[v]->id];
            }
        }

        if (work_mode_ == TE)
        {
            for (int i = 0; i < node_num_; ++i)
            {
                le_.calculate_EH(-1.0*mesh_nodes_[i].y, mesh_nodes_[i].x);
                d1[i] += le_.get_Hy().real();
                d2[i] += le_.get_Hy().imag();
                d3[i] += le_.get_Hz().real();
                d4[i] += le_.get_Hz().imag();
            }

            meshio_.save_data("Hy_real", d1, gctl::NodeData);
            meshio_.save_data("Hy_imag", d2, gctl::NodeData);
            meshio_.save_data("Hz_real", d3, gctl::NodeData);
            meshio_.save_data("Hz_imag", d4, gctl::NodeData);
        }
        else // TM
        {
            for (int i = 0; i < node_num_; ++i)
            {
                le_.calculate_EH(-1.0*mesh_nodes_[i].y, mesh_nodes_[i].x);
                d1[i] += le_.get_Ey().real();
                d2[i] += le_.get_Ey().imag();
                d3[i] += le_.get_Ez().real();
                d4[i] += le_.get_Ez().imag();
            }

            meshio_.save_data("Ey_real", d1, gctl::NodeData);
            meshio_.save_data("Ey_imag", d2, gctl::NodeData);
            meshio_.save_data("Ez_real", d3, gctl::NodeData);
            meshio_.save_data("Ez_imag", d4, gctl::NodeData);
        }
    }

    std::clog << "done\n";
    return;
}

void mt2d::write_observation(std::string filename)
{
    std::clog << "Writing Observation file ... ";
    std::ofstream outfile;
    gctl::open_outfile(outfile, filename, ".obs.txt");
    for (size_t i = 0; i < obs_num_; i++)
    {
        outfile << obs_site_[i].vert->x << " " << obs_site_[i].vert->y << " " << std::setprecision(16)
            << obs_site_[i].app_resist << " " << obs_site_[i].phase << std::endl;
    }
    outfile.close();
    std::clog << "done\n";
    return;
}

double mt2d::get_element_conductivity(int tag_num, bool ref_cndt)
{
    for (int i = 0; i < secondary_model_.size(); ++i) // 数量有限直接搜索即可
    {
        if (secondary_model_[i].tag == tag_num)
        {
            if (ref_cndt) return secondary_model_[i].ref_cndt;
            else return secondary_model_[i].cndt;
        }
    }

    throw gctl::runtime_error("No secondary model found by the tag number: "+std::to_string(tag_num));
    return 0.0;
}

void mt2d::assemble_linear_system()
{
    std::clog << "Assembling kernel matrix ... ";

    std::vector<triplt_cd> kernel_entries;
    kernel_entries.reserve(edge_num_*2);

    g2d_tri_.initiate(2); // 初始化高斯积分阶数

    // 计算线性方程组的左端项
    A.resize(edge_num_, edge_num_);
    A.setZero();

    edge_setup e_set;
    double edge_cndt, aij_p1, aij_p2;
    complex_d tmp_entry;
    gctl::vertex2dc* tmp_vert[3];
    // Loop all elements
    for (int e = 0; e < elem_num_; ++e)
    {
        tmp_vert[0] = mesh_elems_[e].vert[0];
        tmp_vert[1] = mesh_elems_[e].vert[1];
        tmp_vert[2] = mesh_elems_[e].vert[2];

        for (int i = 0; i < 3; i++)
        {
            e_set.uv_id[0] = i;
            if (tmp_vert[i]->id < tmp_vert[(i+1)%3]->id) e_set.uv_type[0] = efem_sf::Normal;
            else e_set.uv_type[0] = efem_sf::Reverse; // determine global edge orientation (from small to big indexing number)

            edge_cndt = get_element_conductivity(elem_markers_[mesh_elems_[e].id]);

            for (int j = 0; j < 3; j++)
            {
                e_set.uv_id[1] = j;
                if (tmp_vert[j]->id < tmp_vert[(j+1)%3]->id) e_set.uv_type[1] = efem_sf::Normal;
                else e_set.uv_type[1] = efem_sf::Reverse;

                aij_p1 = g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Aij_part1, &e_set, this);

                aij_p2 = g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Aij_part2, &e_set, this);

                if (work_mode_ == TE)
                {
                    tmp_entry.real(aij_p1/edge_cndt);
                    tmp_entry.imag(-1.0*omega_*_mu0*aij_p2);
                }
                else
                {
                    tmp_entry.real(aij_p1);
                    tmp_entry.imag(-1.0*omega_*_mu0*edge_cndt*aij_p2);
                }

                kernel_entries.push_back(triplt_cd(elem_edge_ptr_[e][i]->id, elem_edge_ptr_[e][j]->id, tmp_entry));
            }
        }
    }

    A.setFromTriplets(kernel_entries.begin(), kernel_entries.end());
    kernel_entries.clear();

    X.resize(edge_num_);
    X.setZero();

    // 计算线性方程组的右端项 对于一次场方法右端项为0
    B.resize(edge_num_);
    B.setZero();

    // 计算线性方程组的右端项 仅在使用二次场方法时计算
    if (work_mode_ == TM)
    {
        for (int e = 0; e < elem_num_; ++e)
        {
            tmp_vert[0] = mesh_elems_[e].vert[0];
            tmp_vert[1] = mesh_elems_[e].vert[1];
            tmp_vert[2] = mesh_elems_[e].vert[2];

            edge_cndt = get_element_conductivity(elem_markers_[mesh_elems_[e].id]) - get_element_conductivity(elem_markers_[mesh_elems_[e].id], true);

            if (edge_cndt != 0.0) // only calculate elements that have a residual conductivity
            {
                for (int i = 0; i < 3; i++)
                {
                    e_set.uv_id[0] = i;
                    if (tmp_vert[i]->id < tmp_vert[(i+1)%3]->id) e_set.uv_type[0] = efem_sf::Normal;
                    else e_set.uv_type[0] = efem_sf::Reverse;

                    B[elem_edge_ptr_[e][i]->id] += _1i*omega_*_mu0*edge_cndt*g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, 
                        tmp_vert[2]->x, tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Bi_TM, &e_set, this);
                }
            }
        }
    }
    else // work_mode_ == TE
    {
        for (int e = 0; e < elem_num_; ++e)
        {
            tmp_vert[0] = mesh_elems_[e].vert[0];
            tmp_vert[1] = mesh_elems_[e].vert[1];
            tmp_vert[2] = mesh_elems_[e].vert[2];

            edge_cndt = (1.0 - get_element_conductivity(elem_markers_[mesh_elems_[e].id], true)/get_element_conductivity(elem_markers_[mesh_elems_[e].id]));

            if (edge_cndt != 0.0) // only calculate elements that have a residual conductivity
            {
                for (int i = 0; i < 3; i++)
                {
                    e_set.uv_id[0] = i;
                    if (tmp_vert[i]->id < tmp_vert[(i+1)%3]->id) e_set.uv_type[0] = efem_sf::Normal;
                    else e_set.uv_type[0] = efem_sf::Reverse;

                    B[elem_edge_ptr_[e][i]->id] += edge_cndt*g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, 
                        tmp_vert[2]->x, tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Bi_TE, &e_set, this);
                }
            }
        }
    }

    std::clog << "done\n";
    return;
}

void mt2d::apply_boundary_condition()
{
    std::clog << "Applying boundary condition ... ";
    
    for (int i = 0; i < edge_num_; ++i)
    {
        if (mesh_edges_[i].neigh[1] ==  nullptr) // 边界上的顶点
        {
            B.coeffRef(i) = _zero; // 二次场方法计算中假设边界上的二次场已经衰减为0
        }
    }
    
    // 采用划0置1的方法处理第一类边界条件 （稀疏矩阵处理起来很快）

    for (int i = 0; i < edge_num_; ++i) // 外层是列循环
    {
        if (mesh_edges_[i].neigh[1] ==  nullptr)
        {
            for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it) // 只循环非零值 行循环
            {
                if (mesh_edges_[it.row()].neigh[1] != nullptr)  // 非边界上的顶点
                {
                    B[it.row()] -= it.value() * B[i];
                }
            }
        }
    }

    for (int i = 0; i < edge_num_; i++) // 处理行
    {
        for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it)
        {
            if (mesh_edges_[it.row()].neigh[1] == nullptr)
            {
                it.valueRef() = _zero;
            }
        }
    }

    for (int i = 0; i < edge_num_; ++i) // 处理列
    {
        if (mesh_edges_[i].neigh[1] == nullptr)
        {
            for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it)
            {
                it.valueRef() = _zero;
                if (it.row() == it.col()) it.valueRef() = _one; // 置1
            }
        }
    }

    A.makeCompressed();
    // 边界条件处理完成

    std::clog << "done\n";
    return;
}

void mt2d::solve()
{
    std::clog << "Solving the electromagnetic field ... ";
    // Call direct solver to solve the linear system
    Eigen::SparseLU<spmat_cd> solver;
    solver.analyzePattern(A);
    solver.factorize(A);
    X = solver.solve(B);
    std::clog << "done\n";
    return;
}

void mt2d::post_process()
{
    // Sort out MT stations
    mt_station tmp_station;
    for (size_t i = 0; i < node_num_; i++)
    {
        if (node_markers_[i] == station_marker_)
        {
            tmp_station.vert = mesh_nodes_.get(i);
            tmp_station.app_resist = 0.0;
            tmp_station.phase = 0.0;

            obs_site_.push_back(tmp_station);
        }
    }

    obs_num_ = obs_site_.size();

    // Find triangles that include mt stations
    std::vector<mt_station>::iterator stat_iter;
    for (size_t i = 0; i < elem_num_; i++)
    {
        if (elem_markers_[i] != air_marker_) // exclude the air layer to obtain smoother results
        {
            for (size_t j = 0; j < 3; j++)
            {
                if (node_markers_[mesh_elems_[i].vert[j]->id] == station_marker_)
                {
                    for (size_t k = 0; k < obs_num_; k++)
                    {
                        if (obs_site_[k].vert == mesh_elems_[i].vert[j])
                        {
                            obs_site_[k].host_tris.push_back(mesh_elems_.get(i));
                            break;
                        }
                    }
                }
            }
        }
    }

    // Calculate observations
    unsigned int edge_id;
    double total_area, edge_cndt, ux, uy;
    complex_d Ex, Ey, Hx, Hy;
    efem_sf::ordertype_e o_type;
    gctl::vertex2dc* tmp_vert[3];

    if (work_mode_ == TE)
    {
        for (size_t i = 0; i < obs_num_; i++)
        {
            le_.calculate_EH(-1.0*obs_site_[i].vert->y, obs_site_[i].vert->x);

            total_area = 0.0;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                total_area += obs_site_[i].host_tris[j]->area();
            }

            Ex = Hy = _zero;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                edge_cndt = get_element_conductivity(elem_markers_[obs_site_[i].host_tris[j]->id]);

                tmp_vert[0] = obs_site_[i].host_tris[j]->vert[0];
                tmp_vert[1] = obs_site_[i].host_tris[j]->vert[1];
                tmp_vert[2] = obs_site_[i].host_tris[j]->vert[2];

                // 下面其实就是一个平面函数插值过程
                for (size_t k = 0; k < 3; k++)
                {
                    if (tmp_vert[k]->id < tmp_vert[(k+1)%3]->id) o_type = efem_sf::Normal;
                    else o_type = efem_sf::Reverse; // determine global edge orientation (from small to big indexing number)

                    edge_id = elem_edge_ptr_[obs_site_[i].host_tris[j]->id][k]->id;

                    elsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, 
                        k, efem_sf::Value, o_type, ux, uy);

                    Hy += X[edge_id]*ux*(obs_site_[i].host_tris[j]->area()/total_area);

                    elsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, 
                        k, efem_sf::Curl, o_type, ux, uy);

                    Ex += -1.0*(1.0/edge_cndt)*(obs_site_[i].host_tris[j]->area()/total_area)*X[edge_id]*ux;
                }
            }

            // Add the primary field
            Hy += le_.get_Hy();
            Ex += le_.get_Ex();

            obs_site_[i].app_resist = std::norm(Ex/Hy)/(omega_*_mu0);
            obs_site_[i].phase = std::arg(Ex/Hy)*180.0/pi;
        }
    }
    else // work_mode_ == TM
    {
        for (size_t i = 0; i < obs_num_; i++)
        {
            le_.calculate_EH(-1.0*obs_site_[i].vert->y, obs_site_[i].vert->x);

            total_area = 0.0;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                total_area += obs_site_[i].host_tris[j]->area();
            }

            Hx = Ey = _zero;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                tmp_vert[0] = obs_site_[i].host_tris[j]->vert[0];
                tmp_vert[1] = obs_site_[i].host_tris[j]->vert[1];
                tmp_vert[2] = obs_site_[i].host_tris[j]->vert[2];

                // 下面其实就是一个平面函数插值过程
                for (size_t k = 0; k < 3; k++)
                {
                    if (tmp_vert[k]->id < tmp_vert[(k+1)%3]->id) o_type = efem_sf::Normal;
                    else o_type = efem_sf::Reverse; // determine global edge orientation (from small to big indexing number)

                    edge_id = elem_edge_ptr_[obs_site_[i].host_tris[j]->id][k]->id;

                    elsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, 
                        k, efem_sf::Value, o_type, ux, uy);

                    Ey += X[edge_id]*ux*(obs_site_[i].host_tris[j]->area()/total_area);

                    elsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, 
                        k, efem_sf::Curl, o_type, ux, uy);

                    Hx += -1.0*X[edge_id]*ux*(obs_site_[i].host_tris[j]->area()/total_area);
                }
            }
            Hx /= _1i*omega_*_mu0;

            // Add the primary field
            Hx += le_.get_Hx();
            Ey += le_.get_Ey();

            obs_site_[i].app_resist = std::norm(Ey/Hx)/(omega_*_mu0);
            obs_site_[i].phase = std::arg(Ey/Hx)*180.0/pi;
        }
    }

    return;
}

void mt2d::run(int argc, char *argv[])
{
    // Read configuration file
    read_config(argv[1]);
    // Read model file
    read_model(tri_file_);

    if (work_mode_ != Modeling)
    {
        // Prepare the linear system
        assemble_linear_system();
        // Apply boundary condition
        apply_boundary_condition();
        // Solve the linear system
        solve();
        // Calculate observations
        post_process();
        // Save Gmsh model file
        write_model(tri_file_);
        // Write observations
        write_observation(tri_file_);
    }
    else
    {
        // Save Gmsh model file
        write_model(tri_file_, false);
    }
    return;
}

int main(int argc, char *argv[]) try
{
    if (argc < 2)
    {
        std::clog << "Usage: ./mt2d <config_file>\n";
        return -1;
    }

    mt2d instance;
    instance.run(argc, argv);
    return 0;
}
catch (std::exception &e)
{
    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
}
